function [camA xfmAview R TT] = GetCameraS03Matrix( pos, lookat ) 
dbstop if error 
% the matrix that transforms the 
% camera to be at pos, viewing the origin

% The T part of the matrix is just [pos]

camA = zeros(4,4); % we get a 4x4 matrix operating on homog. coordinates

lookat = lookat(:);
T = pos(:);

dir = (lookat-T) / norm(T-lookat); % direction from point to origin

% R^-1 maps dir to [0 0 1]'

phi = atan2( dir(1), dir(2) );
P1 = [cos(phi) -sin(phi) 0; sin(phi) cos(phi) 0; 0 0 1];
dhat = P1*dir;
alpha = pi/2 - atan2(dir(3),dhat(2));
P2 = [1 0 0; 0 cos(alpha) -sin(alpha); 0 sin(alpha) cos(alpha)];

Rinv = P2*P1;
R = [ Rinv zeros(3,1); zeros(1,3), 1];
TT = [ eye(3,3) T; zeros(1,3) 1];
camA = R * TT; % how the camera sees the 'look at' point in its coord system 
xfmAview = TT * R; % transform of the camera itself

% R after TT gives "what is the coordinate of the operated
% upon point in the camera coordinates, where the camera is at T 
% and pointed at the origin



% draw a plot for debug purposes

x_ = [ [0;0;0;1] [1;0;0;1] ];
y_ = [ [0;0;0;1] [0;1;0;1] ];
z_ = [ [0;0;0;1] [0;0;1;1] ];
figure(1); plot3(x_(1,:),x_(2,:),x_(3,:),'r');  hold on
plot3(y_(1,:),y_(2,:),y_(3,:),'g'); plot3(z_(1,:),z_(2,:),z_(3,:),'b');
x_c = camA * x_; y_c = camA * y_; z_c = camA * z_;
figure(1); plot3(x_c(1,:),x_c(2,:),x_c(3,:),'r');  
plot3(y_c(1,:),y_c(2,:),y_c(3,:),'g'); plot3(z_c(1,:),z_c(2,:),z_c(3,:),'b');
hold off

breakhere = 1;